DALS023-机器学习02-条件概率与Loess拟合

发表于 | 分类于
| 字数统计: 3,211 | 阅读时长 ≈ 14

前言

这一部分是《Data Analysis for the life sciences》的第9章机器学习的第2小节，这一部分的主要内容涉及条件概率与Loess回归，这一部分相关的Rmarkdown文档参见作者的Github。

条件概率与期望

预测问题的结果可以分为两类，分别为分类结果(categorical outcomes)和连续型结果(continuous outcomes)。但是，许多算法都可以应用于这两种情况，这是因为条件概率(conditional probabilities)和条件期望(conditional expectations)之间有着一定联系。

对于分类数据，例如二分类结果，如果我们知道给在一组预测因子$X=(X_1,\dots,X_p)^\top$下预测的概率$Y$可能是$k$个结果中的任意一个，这种情况用方程表示就是：

我们就可以优化我们的预测结果。尤其是，对于任意的 $x$ ，我们可能预测出最大概率 $f_k(x)$ 时的 $k$ 。

为了简化我们的描述，我们只考虑二分类变量。我们可以认为 $\mbox{Pr}(Y=1 \mid X=x)$ 是当 $X=x$ 时，第1层的某结果在总体中的比例。考虑到期望是所有 $Y$ 值的均值，在这个案例中，期望的概率就等于：$f(x) \equiv \mbox{E}(Y \mid X=x)=\mbox{Pr}(Y=1 \mid X=x)$。因此，我们在后面的描述中，就会使用期望，因为这种表示更加普遍。

通常来说，期望值是有着比较受欢迎的数学属性，因为它能够缩小预测值 $\hat{Y}$ 和 $Y$ 之间的距离：

预测中的回归问题

我们在前面使用了父子的身高数据来介绍了回归，这里我们使用机器学习来解释一下这个数据分析过程。在我们前面的那个案例中，我们试图通常父母的身高$X$来预测儿子的身高$Y$。这里我们只有一个预测值(predictor)。现在如果我们被问到随机选择一个儿子，这个儿子的身高是多少，我们也许会使用平均值来回答，如下所示：

1
2
3
4
5
6
7
8
library(rafalib)
mypar(1,1)
library(UsingR)
data("father.son")
x=round(father.son$fheight) ##round to nearest inch
y=round(father.son$sheight)
hist(y,breaks=seq(min(y),max(y)))
abline(v=mean(y),col="red",lwd=2)

在这个案例中，我们也可以使用$Y$的分布近似服从正态分布来进行预测，那么也就是说，当我们回答平均值的时候，回答对的概率最大。

现在我们再来考虑一些情况，当我们有了更多的信息后，如何进行预测。例如，我们被告之，这个随机选择的这个儿子的父亲身高是71英寸（高于均值1.25个SD）。那么我们如何预测？

1
2
3
4
mypar(1,2)
plot(x,y,xlab="Father's height in inches",ylab="Son's height in inches",main=paste("correlation =",signif(cor(x,y),2)))
abline(v=c(-0.35,0.35)+71,col="red")
hist(y[x==71],xlab="Heights",nc=8,main="",xlim=range(y))

上图的左侧是儿子身高的散点图，其中红色围起来的部分是给定父亲身高是71英寸这个数据后，对应的儿子的身高。条件分布：儿子的身高数据分布就是被父亲71英寸这个数据限定了。

有了父亲身高是71英寸这个数据后，那么我们来猜儿子的身高时，最好的回答还是期望值(expectation)，但是，我们的数据层(strata)已经发生了改变，也就是说，我们来猜测儿子身高$Y$时，要考虑到限制因素，即$Y=71$。因此，我们可以在这个条件之上，再来计算均值，这个均值就是条件期望(conditional expectation)，因此我们对于任意$x$值来预测时，公式如下所示：

通过微积分我们可以发现，这个结果更加接近于二元正态分布，我们可以换成下面的表示形式：

如果我们使用样本的数据来估计这5个参数，那么我们就会得到回归线：

regression, fig.cap
1
2
3
4
5
6
7
8
mypar(1,2)
plot(x,y,xlab="Father's height in inches",ylab="Son's height in inches",
     main=paste("correlation =",signif(cor(x,y),2)))
abline(v=c(-0.35,0.35)+71,col="red")
fit <- lm(y~x)
abline(fit,col=1)
hist(y[x==71],xlab="Heights",nc=8,main="",xlim=range(y))
abline(v = fit$coef[1] + fit$coef[2]*71, col=1)

在这个特殊的案例中，回归线提供了对 $Y$ 更优化的预测函数，但是通常情况下这种情况并不真实，因为在典型的机器学习问题中，优化后的 $f(x)$ 很少是一条直线。

练习

P379

平滑

相关的Rmarkdown文档见作者的Github。

在所有的数据分析中，平滑(Smoothing)是一个非常强大的工具。当数据的分布的形状未知时，我们可以假定这些数据是平滑(smooth)的，我们就能够估计 $f(x)$ 。平滑的主要思路就是将那些有着类似期望的数据归类，然后计算这些归类后的平均值，或者是对其进行简单的参数模型拟合。我们使用基因表达数据来介绍两个平滑工具。

下面的数据是源于相同RNA样本的基因表达数据。

1
2
3
4
5
6
7
##Following three packages are available from Bioconductor
library(Biobase)
# BiocManager::install("SpikeIn")
library(SpikeIn)
# BiocManager::install("hgu95acdf")
library(hgu95acdf)
data(SpikeIn95)

我们可以使用MA图 ($Y$ = log ratios and $X$ = averages) 来比较这两个重复样本的质量，我们可能通过一种方式对这些数据进行采样，这些采样方式要能平衡 $X$ 中不同层的数据点的数目，如下所示：

echo
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
##Example with two columns
i=10;j=9
##remove the spiked in genes and take random sample
library(GEOquery)
siNames<-colnames(pData(SpikeIn95))
ind <- which(!probeNames(SpikeIn95)%in%siNames)
pms <- pm(SpikeIn95)[ ind ,c(i,j)]
##pick a representative sample for A and order A
Y=log2(pms[,1])-log2(pms[,2])
X=(log2(pms[,1])+log2(pms[,2]))/2
set.seed(4)
ind <- tapply(seq(along=X),round(X*5),function(i)
  if(length(i)>20) return(sample(i,20)) else return(NULL))
ind <- unlist(ind)
X <- X[ind]
Y <- Y[ind]
o <-order(X)
X <- X[o]
Y <- Y[o]

采样结束后，我们来绘制散点图，如下所示：

MAplot, fig.cap
1
2
3
library(rafalib)
mypar()
plot(X,Y)

在上面的MA图中我们可以看到 $Y$ 取决于 $X$。这种取决关系存在着偏差，因为我们发现，它们在重复性方面有着偏差，这就意味着，如果不考虑 $X$， 那么 $Y$ 的均值就是0。现在我们想预测 $f(x)=\mbox{E}(Y \mid X=x)$ ，以便于移除假这种偏差。线性回归无法捕捉到表示上图曲线 $f(x)$ 的曲率(curvature)，如下所示：

MAplot_with_regression_line, fig.cap
1
2
3
4
5
mypar()
plot(X,Y)
fit <- lm(Y~X)
points(X,Y,pch=21,bg=ifelse(Y>fit$fitted,1,3))
abline(fit,col=2,lwd=4,lty=2)

微区间平滑(Bin Smoothing)

对于上述的数据，我们无法使用直线进行拟合，此时我们就要使用平滑处理的思想，也就是说把这些点分散到不同的组里，计算这些不同组的均值，这就是所谓的微区间平滑处理（我没有找到bin smoothing相应的中文翻译，这里说的微区间平滑是我自己造的词）。这种处理的通常思路就是，我们假定这些数据点分为很多微区间(bin)，这样拟合的曲线是足够”平滑“的，那么在这个微区间中的这个曲线就会近似于常数(approximately constant)。如果我们假设这个曲线是常数(constant)的，那么所有位于某个微区间(bin)中的 $Y$ 就有了相同的期望值。例如，在下面的图形中，如果我们将微区间的宽度设为1，我们标记出了在8.6和12.1处的微区间。我们还显示了，在这两个区间中拟合的 $Y$ 的均值，如下所示：

binsmoother, fig.cap
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
mypar()
centers <- seq(min(X),max(X),0.1)
plot(X,Y,col="grey",pch=16)
windowSize <- .5
i <- 25
center<-centers[i]
ind=which(X>center-windowSize & X<center+windowSize)
fit<-mean(Y)
points(X[ind],Y[ind],bg=3,pch=21)
lines(c(min(X[ind]),max(X[ind])),c(fit,fit),col=2,lty=2,lwd=4)
i <- 60
center<-centers[i]
ind=which(X>center-windowSize & X<center+windowSize)
fit<-mean(Y[ind])
points(X[ind],Y[ind],bg=3,pch=21)
lines(c(min(X[ind]),max(X[ind])),c(fit,fit),col=2,lty=2,lwd=4)

我们通过计算每一个点附近的小区间，我们就能对构成的曲线的函数 $f(x) $ 进行估计。下图显示了这个计算过程，即我们从最小的 $x$ 值扩展到最大值过程中的这个估计值，同时我们还展示了最小值与最大值之间的10个值，一共是10张图片，如下所示：

bin_smoothing_demo, fig.cap
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
windowSize<-0.5
smooth<-rep(NA,length(centers))
mypar (4,3)
for(i in seq(along=centers)){
  center<-centers[i]
  ind=which(X>center-windowSize & X<center+windowSize)
  smooth[i]<-mean(Y[ind])
  if(i%%round(length(centers)/12)==1){ ##we show 12
    plot(X,Y,col="grey",pch=16)
    points(X[ind],Y[ind],bg=3,pch=21)
    lines(c(min(X[ind]),max(X[ind])),c(smooth[i],smooth[i]),col=2,lwd=2)
    lines(centers[1:i],smooth[1:i],col="black")
    points(centers[i],smooth[i],col="black",pch=16,cex=1.5)
  }
}

最终的结果就如下所示：

bin_smooth_final, fig.cap
1
2
3
mypar (1,1)
plot(X,Y,col="darkgrey",pch=16)
lines(centers,smooth,col="black",lwd=3)

在R中有许多函数可以进行平滑处理(bin smoothers)，例如ksmooth。但是在实际计算过程中，我们更偏向使用一些稍微复杂的模型来对数据进行拟合。在最后一个案例，也就是最后一张图，这个曲线就不怎么平滑，有一些粗糙。我们后面会介绍loess方法会改善这一点。

Loess

Loess的全称是Local weighted regression，即局部权重回归，它在原理上类似于微区间平滑处理。但Loess与微区间平滑处理的主要区别就在于，在微区间中我是使用直线，还是抛物线进行拟合。Loess计算会增大微区间的数目，但这会使我们的估计更稳定。下图展示了两个微区间的拟合曲线，这两个区间比较宽，使用宽区间主要是因为拟合曲线提供了更多的灵活性，如下所示：

loess, fig.cap
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
centers <- seq(min(X),max(X),0.1)
mypar (1,1)
plot(X,Y,col="darkgrey",pch=16)
windowSize <- 1.25
i <- 25
center<-centers[i]
ind=which(X>center-windowSize & X<center+windowSize)
fit<-lm(Y~X,subset=ind)
points(X[ind],Y[ind],bg=3,pch=21)
a <- min(X[ind]);b <- max(X[ind])
lines(c(a,b),fit$coef[1]+fit$coef[2]*c(a,b),col=2,lty=2,lwd=3)
i <- 60
center<-centers[i]
ind=which(X>center-windowSize & X<center+windowSize)
fit<-lm(Y~X,subset=ind)
points(X[ind],Y[ind],bg=3,pch=21)
a <- min(X[ind]);b <- max(X[ind])
lines(c(a,b),fit$coef[1]+fit$coef[2]*c(a,b),col=2,lty=2,lwd=3)

当我们

现在我们展示一下通过12步处理来对这些数据进行loess拟合，如下所示：

loess_demo, fig.cap
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
library(rafalib)
mypar (4,3)
windowSize<-1.25
smooth<-rep(NA,length(centers))
for(i in seq(along=centers)){
  center<-centers[i]
  ind=which(X>center-windowSize & X<center+windowSize)
  fit<-lm(Y~X,subset=ind)
  smooth[i]<-fit$coef[1]+fit$coef[2]*center
  if(i%%round(length(centers)/12)==1){ ##we show 12
    plot(X,Y,col="grey",pch=16)
    points(X[ind],Y[ind],bg=3,pch=21)
    a <- min(X[ind]);b <- max(X[ind])
    lines(c(a,b),fit$coef[1]+fit$coef[2]*c(a,b),col=2,lwd=2)
    
    lines(centers[1:i],smooth[1:i],col="black")
    points(centers[i],smooth[i],col="black",pch=16,cex=1.5)
  }
}

最终的结果就是生成一条更加平滑的拟合曲线用于估计我们的局部参数，如下所示：

loess_final, fig.cap
1
2
3
mypar (1,1)
plot(X,Y,col="darkgrey",pch=16)
lines(centers,smooth,col="black",lwd=3)

R中的函数loess()可以进行上述的分析，如下所示：

loess2, fig.cap
1
2
3
4
5
6
fit <- loess(Y~X, degree=1, span=1/3)
newx <- seq(min(X),max(X),len=100) 
smooth <- predict(fit,newdata=data.frame(X=newx))
mypar ()
plot(X,Y,col="darkgrey",pch=16)
lines(newx,smooth,col="black",lwd=3)

loess()函数与其它的曲线平滑处理函数(bin smoother)还有三点重要的区别：

第一，loess()函数会保持局部拟合中的点数不变，而非保持区间的数目不变。局部拟合中的点数是通过span参数决定的，这是一个比值。例如，如果 N 是数据点的数目，那么 span = 0.5 则表示针对一个确定的 x ， loess()将会使用 0.5*N 个最接近的点对 x 进行拟合。

第二，使用参数模型来对 $f(x)$ 进行拟合时，loess()会使用加权最小平方和来进行计算，那些权重比较大的点就是那些更接近 x 的点。

第三，loess()函数会选择更稳健(robustly)的局部模型来进行拟合。这是一种迭代算法，其中在一次迭代中拟合了模型之后，对于那些检测到的异常值就会在下一次迭代中降低其权重。如果要想使用这个参数，可以设定family = "symmetric"。

练习

P389